module mardiamodadj
implicit none
contains
subroutine mardiajbfadj(m,T,n,y,theta,testmardia,testjb)
use linear_operators
use matrixfuns
use datahadler
use mardiamod
integer, intent(in):: m,T,n
integer, parameter:: p=11,pjb=2*3
double precision, intent(in):: theta(m),y(T,n)
double precision, intent(out):: testmardia(3),testjb(3)
integer i,j,ii,irank
double precision, parameter:: c9=.999d0
double precision lamb(T),delta(n),c(n),phi0(n),phi1(n),phi2(n),gamma(n),&
                &sigmat(n,n),gammat(n,n),alpha1,alpha2,vc(n,1),&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n)
double precision ddelta(m,n),dc(m,n),dgamma(m,n),ent(n,n*n)&
               &,dmut(m,n),dalpha1(m),dalpha2(m),alpha0,dalpha0(m)&
			   &,dlambda(m,1),ident(n,n),wttm,fttm&
			   &,dfttm(m),dwttm(m),matwtt(m),metn(n),dphi0(m,n),dphi1(m,n),dphi2(m,n)
double precision dvsig(m,n*n),ckc(m,n*n),ceye(m,n*n),matB(p,m),matR(p-1,n),rsigmat(n,n)&
&,irsigmat(n,n),tol,dmach,matI(m,m),moments(p),epstar(n),wmat(p,p)&
&,momjb(pjb),matBjb(pjb,m),mkujb(n,n*n),mskjb(n,n),wmatjb(pjb,pjb)


ent=matent(n)
ident=eye(n)
delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall

vc(:,1)=c

forall(i=1:m,j=1:n)
	ddelta(i,j)=0.0d0
	dc(i,j)=0.0d0
	dgamma(i,j)=0.0d0
	dphi0(i,j)=0.0d0
	dphi1(i,j)=0.0d0
	dphi2(i,j)=0.0d0
	dmut(i,j)=0.0d0
end forall
forall(i=1:n)
	ddelta(i,i)=1.0d0
	dphi0(2*n+1,i)=1.0d0
	dphi1(2*n+4,i)=c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4))
	dphi2(2*n+4,i)=-(c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4)))*(dsin(theta(2*n+5))**2.0d0)
	dphi2(2*n+5,i)=(c9-phi1(i))*2.0d0*dsin(theta(2*n+5))*dcos(theta(2*n+5))
end forall
forall(i=2:n)
	dc(n+i-1,i)=1.0d0
end forall

forall(i=1:m)
	dalpha0(i)=0.0d0
	dalpha1(i)=0.0d0
	dlambda(i,1)=0.0d0
end forall
dalpha2=dalpha1
dalpha0(2*n)=1.0d0
dalpha1(2*n+2)=c9*2.0d0*dsin(theta(2*n+2))*dcos(theta(2*n+2))
dalpha2(2*n+2)=-dalpha1(2*n+2)*(dsin(theta(2*n+3))**2.0d0)
dalpha2(2*n+3)=(c9-alpha1)*2.0d0*dsin(theta(2*n+3))*dcos(theta(2*n+3))


forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:n)
	dmut(i,i)=1.0d0
end forall

gamma=phi0/(1.0d0-phi1-phi2)

gammat=diag(gamma)
igammat=diag(1.0d0/gamma)
forall (i=1:m,j=1:n)
	dgamma(i,j)=((1.0d0-phi1(j)-phi2(j))*dphi0(i,j)-phi0(j)*(-dphi1(i,j)-dphi2(i,j)))&
	           &/(1.0d0-phi1(j)-phi2(j))**2.0d0
end forall
dlambda(:,1)=((1.0d0-alpha1-alpha2)*dalpha0-alpha0*(-dalpha1-dalpha2))/(1.0d0-alpha1-alpha2)**2.0d0

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)
sigmat=lamb(1)*cc+gammat
isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(1))))


ckc=.t.kron(vc,vc)
ceye=matmul(dc,kron(.t.vc,dble(eye(n)))+kron(dble(eye(n)),.t.vc))
forall(i=1:p-1,j=1:n)
	matR(i,j)=0.0d0
end forall
matR(1,1)=3.0d0
matR(2,2)=3.0d0
matR(3,3)=3.0d0
matR(4,2)=1.0d0
matR(5,3)=1.0d0
matR(6,1)=1.0d0
matR(7,3)=1.0d0
matR(8,1)=1.0d0
matR(9,2)=1.0d0

tol = 100.0*dmach(4)

call dchfac(n,sigmat,n,tol,irank,rsigmat,n)
irsigmat=.i.(.t.rsigmat)

dvsig=lamb(1)*ceye+matmul(dgamma,ent)

matB(1,:)=2.0d0*(n+2.0d0)*&
         &(2.0d0*lamb(1)*matmul(dc,matmul(isigmat,c))+matmul(dgamma,vecd(isigmat))&
		 &+dot_product(c,matmul(isigmat,c))*dlambda(:,1))
matB(2:p,:)=matmul(matR,irsigmat.xt.dmut)
matI=matmul(dmut,isigmat.xt.dmut)+.5d0*matmul(dvsig,kron(isigmat,isigmat).xt.dvsig)

epstar=matmul(irsigmat,epst(1,:))
moments(1)=(dot_product(epstar,epstar)**2.0d0)-n*(n+2.0d0)
moments(2:4)=epstar**3.0d0
moments(5:6)=(epstar(1)**2.0d0)*epstar(2:3)
moments(7)=(epstar(2)**2.0d0)*epstar(1)
moments(8)=(epstar(2)**2.0d0)*epstar(3)
moments(9:10)=(epstar(3)**2.0d0)*epstar(1:2)
moments(11)=epstar(1)*epstar(2)*epstar(3)

!!!!
forall(i=1:n,j=1:n*n)
	mkujb(i,j)=0.0d0
end forall

mkujb(1,1)=12.0d0
mkujb(2,5)=12.0d0
mkujb(3,9)=12.0d0
mskjb=3.0d0*eye(n)

!momjb(1:3)=(epstar**4.0d0)-3.0d0
!momjb(4:6)=(epstar**3.0d0)-3.0d0*epstar
!matBjb(1:n,:)=.5d0*matmul(mkujb,kron(irsigmat,irsigmat).xt.dvsig)
!forall(i=n+1:2*n,j=1:m)
!matBjb(i,j)=0.0d0
!end forall

momjb(1:3)=(epstar**4.0d0)-3.0d0
momjb(4:6)=epstar**3.0d0
matBjb(1:n,:)=.5d0*matmul(mkujb,kron(irsigmat,irsigmat).xt.dvsig)
matBjb(n+1:2*n,:)=matmul(mskjb,irsigmat.xt.dmut)

!t>1
do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0+alpha1*((fttm**2.0d0)+wttm)+alpha2*lamb(i-1)

	forall(ii=1:n)
	    metn(ii)=dot_product(igammat(ii,:),c)*dot_product(igammat(:,ii),epst(i-1,:))
	end forall
	matwtt=matmul(matmul(dgamma,ent),vec2((igammat.x.cc).x.igammat))

	dwttm=-(wttm*wttm)*(-(lamb(i-1)**(-2.0))*dlambda(:,1)+2.0d0*matmul(dc,matmul(igammat,c))&
	      &-matwtt)
	
	dfttm=dwttm*dot_product(c,matmul(igammat,epst(i-1,:)))&
	     &+wttm*matmul(dc,matmul(igammat,epst(i-1,:)))&
		 &-wttm*matmul(dgamma,metn)-wttm*matmul(dmut,matmul(igammat,c))


	forall (ii=1:m,j=1:n)
		dgamma(ii,j)=dphi0(ii,j)+dphi1(ii,j)*(((epst(i-1,j)-c(j)*fttm)**2.0d0)&
		                                                                   &+(c(j)*c(j))*wttm)&
		           &+dphi2(ii,j)*gamma(j)&
				   &+phi1(j)*2.0d0*(epst(i-1,j)-c(j)*fttm)*&
				   &(-dmut(ii,j)-dc(ii,j)*fttm-c(j)*dfttm(ii))&
				   &+phi1(j)*2.0d0*c(j)*dc(ii,j)*wttm+phi1(j)*c(j)*c(j)*dwttm(ii)+phi2(j)*dgamma(ii,j)
	end forall


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)
	igammat=diag(1.0d0/gamma)


	sigmat=lamb(i)*cc+gammat
	isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(i))))

	dlambda(:,1)=dalpha0+dalpha1*((fttm*fttm)+wttm)+dalpha2*(lamb(i-1))&
	      &+alpha1*(2.0d0*fttm*dfttm+dwttm)+alpha2*dlambda(:,1)

	call dchfac(n,sigmat,n,tol,irank,rsigmat,n)
	irsigmat=.i.(.t.rsigmat)

	dvsig=lamb(i)*ceye+matmul(dgamma,ent)+matmul(dlambda,ckc)
	matB(1,:)=matB(1,:)+ ( 2.0d0*(n+2.0d0)*&
         &(2.0d0*lamb(i)*matmul(dc,matmul(isigmat,c))+matmul(dgamma,vecd(isigmat))&
		 &+dot_product(c,matmul(isigmat,c))*dlambda(:,1)) )
	matB(2:p,:)=matB(2:p,:)+matmul(matR,irsigmat.xt.dmut)
	matI=matI+(matmul(dmut,isigmat.xt.dmut)+.5d0*matmul(dvsig,kron(isigmat,isigmat).xt.dvsig))

	epstar=matmul(irsigmat,epst(i,:))
	moments(1)=moments(1)+((dot_product(epstar,epstar)**2.0d0)-n*(n+2.0d0))
	moments(2:4)=moments(2:4)+(epstar**3.0d0)
	moments(5:6)=moments(5:6)+((epstar(1)**2.0d0)*epstar(2:3))
	moments(7)=moments(7)+((epstar(2)**2.0d0)*epstar(1))
	moments(8)=moments(8)+((epstar(2)**2.0d0)*epstar(3))
	moments(9:10)=moments(9:10)+((epstar(3)**2.0d0)*epstar(1:2))
	moments(11)=moments(11)+(epstar(1)*epstar(2)*epstar(3))

	!!!
	!momjb(1:3)=momjb(1:3)+((epstar**4.0d0)-3.0d0)
	!momjb(4:6)=momjb(4:6)+(epstar**3.0d0)-(3.0d0*epstar)
	!matBjb(1:n,:)=matBjb(1:n,:)+(.5d0*matmul(mkujb,kron(irsigmat,irsigmat).xt.dvsig))


	momjb(1:3)=momjb(1:3)+((epstar**4.0d0)-3.0d0)
	momjb(4:6)=momjb(4:6)+(epstar**3.0d0)
	matBjb(1:n,:)=matBjb(1:n,:)+(.5d0*matmul(mkujb,kron(irsigmat,irsigmat).xt.dvsig))
	matBjb(n+1:2*n,:)=matBjb(n+1:2*n,:)+(matmul(mskjb,irsigmat.xt.dmut))

end do
matB=matB/T
matI=matI/T
matBjb=matBjb/T
call var0(wmat)
!!
forall(i=1:pjb,j=1:pjb)
	wmatjb(i,j)=0.0d0
end forall
wmatjb(1,1)=96.0d0
wmatjb(2,2)=96.0d0
wmatjb(3,3)=96.0d0

wmatjb(4,4)=15.0d0
wmatjb(5,5)=15.0d0
wmatjb(6,6)=15.0d0

!!
wmat=wmat-((matB.x.(.i.matI)).xt.matB)
wmatjb=wmatjb-((matBjb.x.(.i.matI)).xt.matBjb)


testmardia(1)=(1.0d0/T)*moments(1)*moments(1)/wmat(1,1)
testmardia(2)=(1.0d0/T)*dot_product(moments(2:p),matmul(.i.wmat(2:p,2:p),moments(2:p)))
testmardia(3)=(1.0d0/T)*dot_product(moments,matmul(.i.wmat,moments))

testjb(1)=(1.0d0/T)*dot_product(momjb(1:n),matmul(.i.wmatjb(1:n,1:n),momjb(1:n)))
testjb(2)=(1.0d0/T)*dot_product(momjb(n+1:2*n),matmul(.i.wmatjb(n+1:2*n,n+1:2*n),momjb(n+1:2*n)))
testjb(3)=(1.0d0/T)*dot_product(momjb,matmul(.i.wmatjb,momjb))

!call exportvec(vec2(matI),"matIe.txt",1)
!call exportvec(vec2(matB.x.(.i.matI)),"matBe.txt",1)
!call exportvec(vec2(wmat),"wmat.txt",1)

end subroutine mardiajbfadj


subroutine mardiajbfadj2(m,T,n,y,theta,testmardia,testjb)
!Moments orthogonalised with respect to theta
use linear_operators
use matrixfuns
use datahadler
use mardiamod
integer, intent(in):: m,T,n
integer, parameter:: p=11,pjb=2*3
double precision, intent(in):: theta(m),y(T,n)
double precision, intent(out):: testmardia(3),testjb(3)
integer i,j,ii,irank
double precision, parameter:: c9=.999d0
double precision lamb(T),delta(n),c(n),phi0(n),phi1(n),phi2(n),gamma(n),&
                &sigmat(n,n),gammat(n,n),alpha1,alpha2,vc(n,1),&
				&cc(n,n),epst(1:T,1:n),isigmat(n,n),igammat(n,n)
double precision ddelta(m,n),dc(m,n),dgamma(m,n),ent(n,n*n)&
               &,dmut(m,n),dalpha1(m),dalpha2(m),alpha0,dalpha0(m)&
			   &,dlambda(m,1),ident(n,n),wttm,fttm&
			   &,dfttm(m),dwttm(m),matwtt(m),metn(n),dphi0(m,n),dphi1(m,n),dphi2(m,n)
double precision dvsig(m,n*n),ckc(m,n*n),ceye(m,n*n),matB(p,m),matR(p-1,n),rsigmat(n,n)&
&,irsigmat(n,n),tol,dmach,matI(m,m),moments(p),epstar(n),wmat(p,p)&
&,momjb(pjb),matBjb(pjb,m),mkujb(n,n*n),mskjb(n,n),wmatjb(pjb,pjb)&
&,grad(m),siget(n),sigeet(n)


ent=matent(n)
ident=eye(n)
delta=theta(1:n)
c(1)=1.0d0
c(2:n)=theta(n+1:2*n-1)
alpha0=theta(2*n)
phi0=vassign(n,theta(2*n+1))
alpha1=c9*(dsin(theta(2*n+2))**2.0d0)
alpha2=(c9-alpha1)*(dsin(theta(2*n+3))**2.0d0)

forall(i=1:n)
phi1(i)=c9*(dsin(theta(2*n+4))**2.0d0)
end forall
forall(i=1:n)
phi2(i)=(c9-phi1(i))*(dsin(theta(2*n+5))**2.0d0)
end forall

vc(:,1)=c

forall(i=1:m,j=1:n)
	ddelta(i,j)=0.0d0
	dc(i,j)=0.0d0
	dgamma(i,j)=0.0d0
	dphi0(i,j)=0.0d0
	dphi1(i,j)=0.0d0
	dphi2(i,j)=0.0d0
	dmut(i,j)=0.0d0
end forall
forall(i=1:n)
	ddelta(i,i)=1.0d0
	dphi0(2*n+1,i)=1.0d0
	dphi1(2*n+4,i)=c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4))
	dphi2(2*n+4,i)=-(c9*2.0d0*dsin(theta(2*n+4))*dcos(theta(2*n+4)))*(dsin(theta(2*n+5))**2.0d0)
	dphi2(2*n+5,i)=(c9-phi1(i))*2.0d0*dsin(theta(2*n+5))*dcos(theta(2*n+5))
end forall
forall(i=2:n)
	dc(n+i-1,i)=1.0d0
end forall

forall(i=1:m)
	dalpha0(i)=0.0d0
	dalpha1(i)=0.0d0
	dlambda(i,1)=0.0d0
end forall
dalpha2=dalpha1
dalpha0(2*n)=1.0d0
dalpha1(2*n+2)=c9*2.0d0*dsin(theta(2*n+2))*dcos(theta(2*n+2))
dalpha2(2*n+2)=-dalpha1(2*n+2)*(dsin(theta(2*n+3))**2.0d0)
dalpha2(2*n+3)=(c9-alpha1)*2.0d0*dsin(theta(2*n+3))*dcos(theta(2*n+3))


forall(i=1:T,j=1:n)
	epst(i,j)=y(i,j)-delta(j)
end forall


forall(i=1:n,j=1:n)
	cc(i,j)=c(i)*c(j)
end forall

forall(i=1:n)
	dmut(i,i)=1.0d0
end forall

gamma=phi0/(1.0d0-phi1-phi2)

gammat=diag(gamma)
igammat=diag(1.0d0/gamma)
forall (i=1:m,j=1:n)
	dgamma(i,j)=((1.0d0-phi1(j)-phi2(j))*dphi0(i,j)-phi0(j)*(-dphi1(i,j)-dphi2(i,j)))&
	           &/(1.0d0-phi1(j)-phi2(j))**2.0d0
end forall
dlambda(:,1)=((1.0d0-alpha1-alpha2)*dalpha0-alpha0*(-dalpha1-dalpha2))/(1.0d0-alpha1-alpha2)**2.0d0

lamb(1)=alpha0/(1.0d0-alpha1-alpha2)
sigmat=lamb(1)*cc+gammat
isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(1))))


ckc=.t.kron(vc,vc)
ceye=matmul(dc,kron(.t.vc,dble(eye(n)))+kron(dble(eye(n)),.t.vc))
forall(i=1:p-1,j=1:n)
	matR(i,j)=0.0d0
end forall
matR(1,1)=3.0d0
matR(2,2)=3.0d0
matR(3,3)=3.0d0
matR(4,2)=1.0d0
matR(5,3)=1.0d0
matR(6,1)=1.0d0
matR(7,3)=1.0d0
matR(8,1)=1.0d0
matR(9,2)=1.0d0

tol = 100.0*dmach(4)

call dchfac(n,sigmat,n,tol,irank,rsigmat,n)
irsigmat=.i.(.t.rsigmat)

dvsig=lamb(1)*ceye+matmul(dgamma,ent)

matB(1,:)=2.0d0*(n+2.0d0)*&
         &(2.0d0*lamb(1)*matmul(dc,matmul(isigmat,c))+matmul(dgamma,vecd(isigmat))&
		 &+dot_product(c,matmul(isigmat,c))*dlambda(:,1))
matB(2:p,:)=matmul(matR,irsigmat.xt.dmut)
matI=matmul(dmut,isigmat.xt.dmut)+.5d0*matmul(dvsig,kron(isigmat,isigmat).xt.dvsig)

epstar=matmul(irsigmat,epst(1,:))
moments(1)=(dot_product(epstar,epstar)**2.0d0)-n*(n+2.0d0)
moments(2:4)=epstar**3.0d0
moments(5:6)=(epstar(1)**2.0d0)*epstar(2:3)
moments(7)=(epstar(2)**2.0d0)*epstar(1)
moments(8)=(epstar(2)**2.0d0)*epstar(3)
moments(9:10)=(epstar(3)**2.0d0)*epstar(1:2)
moments(11)=epstar(1)*epstar(2)*epstar(3)

!!!!
forall(i=1:n,j=1:n*n)
	mkujb(i,j)=0.0d0
end forall

mkujb(1,1)=12.0d0
mkujb(2,5)=12.0d0
mkujb(3,9)=12.0d0
mskjb=3.0d0*eye(n)


momjb(1:3)=(epstar**4.0d0)-3.0d0
momjb(4:6)=epstar**3.0d0
matBjb(1:n,:)=.5d0*matmul(mkujb,kron(irsigmat,irsigmat).xt.dvsig)
matBjb(n+1:2*n,:)=matmul(mskjb,irsigmat.xt.dmut)

forall(i=1:n)
	sigeet(i)=(dot_product(isigmat(i,:),epst(1,:))*&
	          & dot_product(isigmat(:,i),epst(1,:)))-isigmat(i,i)
end forall
siget=matmul(isigmat,epst(1,:))

grad=matmul(dmut,siget)+.5d0*matmul(dgamma,sigeet) &
&+.5d0*dlambda(:,1)*((dot_product(c,siget)**2.0)-dot_product(c,matmul(isigmat,c)))&
&+.5d0*lamb(1)*(dot_product(c,siget)*2.0d0*matmul(dc,siget)&
&-2.0*matmul(matmul(dc,isigmat),c))
!t>1
do i=2,T,1
	wttm=lamb(i-1)/(1.0d0+lamb(i-1)*dot_product(c,c/gamma))
	fttm=wttm*dot_product((c/gamma),epst(i-1,:))
	lamb(i)=alpha0+alpha1*((fttm**2.0d0)+wttm)+alpha2*lamb(i-1)

	forall(ii=1:n)
	    metn(ii)=dot_product(igammat(ii,:),c)*dot_product(igammat(:,ii),epst(i-1,:))
	end forall
	matwtt=matmul(matmul(dgamma,ent),vec2((igammat.x.cc).x.igammat))

	dwttm=-(wttm*wttm)*(-(lamb(i-1)**(-2.0))*dlambda(:,1)+2.0d0*matmul(dc,matmul(igammat,c))&
	      &-matwtt)
	
	dfttm=dwttm*dot_product(c,matmul(igammat,epst(i-1,:)))&
	     &+wttm*matmul(dc,matmul(igammat,epst(i-1,:)))&
		 &-wttm*matmul(dgamma,metn)-wttm*matmul(dmut,matmul(igammat,c))


	forall (ii=1:m,j=1:n)
		dgamma(ii,j)=dphi0(ii,j)+dphi1(ii,j)*(((epst(i-1,j)-c(j)*fttm)**2.0d0)&
		                                                                   &+(c(j)*c(j))*wttm)&
		           &+dphi2(ii,j)*gamma(j)&
				   &+phi1(j)*2.0d0*(epst(i-1,j)-c(j)*fttm)*&
				   &(-dmut(ii,j)-dc(ii,j)*fttm-c(j)*dfttm(ii))&
				   &+phi1(j)*2.0d0*c(j)*dc(ii,j)*wttm+phi1(j)*c(j)*c(j)*dwttm(ii)+phi2(j)*dgamma(ii,j)
	end forall


	gamma=phi0+phi1*(((epst(i-1,:)-c*fttm)**2.0d0)+(c*c)*wttm)+phi2*gamma

	gammat=diag(gamma)
	igammat=diag(1.0d0/gamma)


	sigmat=lamb(i)*cc+gammat
	isigmat=igammat-(((igammat.x.cc).x.igammat)&
	      &/(dot_product(c,matmul(igammat,c))+(1.0d0/lamb(i))))

	dlambda(:,1)=dalpha0+dalpha1*((fttm*fttm)+wttm)+dalpha2*(lamb(i-1))&
	      &+alpha1*(2.0d0*fttm*dfttm+dwttm)+alpha2*dlambda(:,1)

	call dchfac(n,sigmat,n,tol,irank,rsigmat,n)
	irsigmat=.i.(.t.rsigmat)

	dvsig=lamb(i)*ceye+matmul(dgamma,ent)+matmul(dlambda,ckc)
	matB(1,:)=matB(1,:)+ ( 2.0d0*(n+2.0d0)*&
         &(2.0d0*lamb(i)*matmul(dc,matmul(isigmat,c))+matmul(dgamma,vecd(isigmat))&
		 &+dot_product(c,matmul(isigmat,c))*dlambda(:,1)) )
	matB(2:p,:)=matB(2:p,:)+matmul(matR,irsigmat.xt.dmut)
	matI=matI+(matmul(dmut,isigmat.xt.dmut)+.5d0*matmul(dvsig,kron(isigmat,isigmat).xt.dvsig))

	epstar=matmul(irsigmat,epst(i,:))
	moments(1)=moments(1)+((dot_product(epstar,epstar)**2.0d0)-n*(n+2.0d0))
	moments(2:4)=moments(2:4)+(epstar**3.0d0)
	moments(5:6)=moments(5:6)+((epstar(1)**2.0d0)*epstar(2:3))
	moments(7)=moments(7)+((epstar(2)**2.0d0)*epstar(1))
	moments(8)=moments(8)+((epstar(2)**2.0d0)*epstar(3))
	moments(9:10)=moments(9:10)+((epstar(3)**2.0d0)*epstar(1:2))
	moments(11)=moments(11)+(epstar(1)*epstar(2)*epstar(3))

	!!!

	momjb(1:3)=momjb(1:3)+((epstar**4.0d0)-3.0d0)
	momjb(4:6)=momjb(4:6)+(epstar**3.0d0)
	matBjb(1:n,:)=matBjb(1:n,:)+(.5d0*matmul(mkujb,kron(irsigmat,irsigmat).xt.dvsig))
	matBjb(n+1:2*n,:)=matBjb(n+1:2*n,:)+(matmul(mskjb,irsigmat.xt.dmut))


	forall(ii=1:n)
		sigeet(ii)=(dot_product(isigmat(ii,:),epst(i,:))*&
		          &   dot_product(isigmat(:,ii),epst(i,:)))-isigmat(ii,ii)
	end forall
	siget=matmul(isigmat,epst(i,:))
	grad=grad+matmul(dmut,siget)+.5d0*matmul(dgamma,sigeet)&
		&+.5d0*dlambda(:,1)*((dot_product(c,siget)**2.0)-dot_product(c,matmul(isigmat,c)))&
		&+.5d0*lamb(i)*(dot_product(c,siget)*2.0d0*matmul(dc,siget)&
		&-2.0*matmul(matmul(dc,isigmat),c))

end do
matB=matB/T
matI=matI/T
matBjb=matBjb/T
call var0(wmat)
!!
forall(i=1:pjb,j=1:pjb)
	wmatjb(i,j)=0.0d0
end forall
wmatjb(1,1)=96.0d0
wmatjb(2,2)=96.0d0
wmatjb(3,3)=96.0d0

wmatjb(4,4)=15.0d0
wmatjb(5,5)=15.0d0
wmatjb(6,6)=15.0d0

!!
wmat=wmat-((matB.x.(.i.matI)).xt.matB)
wmatjb=wmatjb-((matBjb.x.(.i.matI)).xt.matBjb)

!!

moments=moments-matmul(matmul(matB,.i.matI),grad)
momjb=momjb-matmul(matmul(matBjb,.i.matI),grad)

testmardia(1)=(1.0d0/T)*moments(1)*moments(1)/wmat(1,1)
testmardia(2)=(1.0d0/T)*dot_product(moments(2:p),matmul(.i.wmat(2:p,2:p),moments(2:p)))
testmardia(3)=(1.0d0/T)*dot_product(moments,matmul(.i.wmat,moments))

testjb(1)=(1.0d0/T)*dot_product(momjb(1:n),matmul(.i.wmatjb(1:n,1:n),momjb(1:n)))
testjb(2)=(1.0d0/T)*dot_product(momjb(n+1:2*n),matmul(.i.wmatjb(n+1:2*n,n+1:2*n),momjb(n+1:2*n)))
testjb(3)=(1.0d0/T)*dot_product(momjb,matmul(.i.wmatjb,momjb))
end subroutine mardiajbfadj2

!!
end module mardiamodadj